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A two-dimensional Ginzburg-Landau theory of weak links in a p-wave superfluid is presented. 
First we consider the symmetry properties of the energy functionals, and their relation to the 
conserved supercurrents which play an essential role in the weak link problem. In numerical studies, 
we use the A and B phases of superfluid 3 He. The phases on the two sides of the weak link can be 
chosen separately, and very general soft degrees of freedom may be imposed as boundary conditions. 
We study all four inequivalent combinations of A and B which are possible for a hole in a planar 
wall, including weak links with a pinned A-B interface. In all cases, some illustrative current-phase 
relations (CPR's) are calculated and the critical currents are mapped. Phase diagrams covering 
the relevant phase space in zero magnetic field are constructed. The numerical methods are also 
described in some detail. 
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I. INTRODUCTION 

Recent experimental observations jl], 0] of the Joseph- 
son effect in weak links of superfluid 'Tie left theorists 
with some interesting problems. One of these was re- 
lated to the interpretation of the u ir states" , local minima 
of the Josephson energy at phase differences other than 
(mod 2tt). Inspired by these findings, some theoreti- 
cal work was done.|| EL M, ||, [?], || By now there exists 
a reasonably good understanding of how such n states 
could be obtained in small "pinhole" contacts, 0, 0, J| in 
arrays of apertures via the anisotextural effect, ffO or 
in single large apertures in the Ginzburg-Landau (GL) 
regime.]^, ^| However, large physically relevant parts of 
the parameter space remain unexplored. In this paper 
we study the GL regime in more detail, and attempt 
to bridge the gap between the pinhole and the large- 
aperture limits. 

In Ref. | we studied a three-dimensional (3D) circu- 
lar aperture between two bulk volumes of B phase 3 He 
by solving the GL equations in a full 3D lattice. With 
this calculation we were able to demonstrate that a local 
minimum of energy exists at phase differences close to 
7T, which is associated with a separate il n branch" in the 
current-phase relation. This solution was never clearly 
found to be a global minimum of energy, and no jumps 
to it from the strongly hysteretic "0 branches" could be 
observed. Therefore it appeared to be experimentally 
inaccessible. This calculation was very restricted, how- 
ever, since the bulk order parameters were fixed to be 
equal on the two sides. In the present calculation we are 
able to impose much more general boundary conditions 
on the soft degrees of freedom of the order parameters, 
and to use the A phase in addition to the B phase. On 
the other hand, we have not attempted to carry out the 
calculations in 3D, but rather we use a two-dimensional 
(2D) model of a long slit-shaped weak link in a planar 
wall.pf Nevertheless, as our experience obtained with the 



previous calculation of Ref. |^ suggests, we may expect 
very similar effects to exist in all weak links, regardless 
of their shape. In the A phase this is not quite so, as will 
be discussed below, because the critical current can even 
vanish in some very symmetrical situations. 

Since this article is to appear in the proceedings of 
a winter school, our approach is in many ways tutorial. 
Furthermore, since the involved work is largely numer- 
ical, special attention is paid to explaining the compu- 
tational methods. Throughout, we shall only deal with 
thermodynamic equilibrium, and the numerics are thus 
mostly related to a minimization of the GL free-energy 
functional on a square lattice. Our relaxation methods 
are actually quite standard ones, but have never been 
presented in detail in the specific context of 3 He. 

In Section 2 we discuss the Josephson coupling on a for- 
mal level, and remind how non-sinusoidal energy-phase 
relations can be obtained via a perturbative approach 
in 3 He and unconventional superconductors alike. From 
Section 3 onward we turn to the case of 3 He, and first re- 
view some basic issues and notations. We relate the sym- 
metries of the free energy functional to the conserved su- 
percurrents, and review the GL theory briefly. In partic- 
ular, we consider spin supercurrents on an equal footing 
with the mass supercurrents everywhere. Section 4 in- 
troduces our 2D weak link model, with the divisions into 
numerical and asymptotic regions, and discusses physi- 
cal ways to control the boundary conditions. Finally, an 
analysis of the four different phase combinations is given 
in Section 5, and some examples of current-phase rela- 
tions (CPR's) are presented. All critical currents and 
phase diagrams are summarized in Section 6. In Sec- 
tion 7 we end with some conclusions and discussion. The 
asymptotic solutions far from the weak link are devel- 
oped in Appendix A, and an analysis of the numerical 
method is presented in Appendix B. 
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II. DC JOSEPHSON EFFECT REVISITED 

The superfluid state of a BCS super- 
fluid/superconductor is described by a gap matrix, 
or a pairing potential. This 2x2 matrix has the form 

A , J g(k) = [A s ia y + A • oia y ] a p. (1) 

Here A s and A are the singlet and triplet pairing ampli- 
tudes, respectively, and their argument k parametrizes 
points on the Fermi surface. From the Pauli principle 
A Q; 3(k) = — A / 3 Q (— k) it follows that A s and A have 
even parity (s, d,g, . . .) and odd parity (p, f,h,.. .) in k, 
respectively. The gap matrix also serves as the order 
parameter of the superfluid. 

The dc (or equilibrium) Josephson effect results from a 
coupling of two superfluids (left, L, and right, R) through 
a weak link. Due to this coupling, the free energy of the 
system is generally changed by an amount Fj which we 
call the Josephson coupling energy. We assume Fj to 
be a function of the bulk values A LR of the order pa- 
rameters. From the order parameters we separate out 
phase factors by writing A LR = Aq ' R exp(icf> L ' R ). Be- 
cause of global gauge invariancc Fj should independent 
of the global phase. The dependence of Fj on the phase 
difference A<fi = <j> R — 4> L is called the energy-phase rela- 
tion (EPR). Since the phase factors are 2ir periodic, so is 
the EPR: Fj{A(j> + 2tt) = Fj(A<f>). As is well known, the 
derivative of EPR with respect to A<fi gives the current- 
phase relation (CPR) for mass supercurrent, J s (A^>).jl(J 
In the triplet case there may also exist a spin current, if 
A L |fA fl . 

In non-magnetic situations it is also reasonable to as- 
sume a symmetry of Fj(A(/)) under time-reversal (TR). 
If A L,R are both "TR invariant", meaning that Ay' R can 
be chosen real so that TR only complex conjugates the 
phase factor, then Fj(-A<j)) = Fj(A0).pjj Similarly, 
the CPR then satisfies J s (— A0) = — J s [A(f>). Let us 
apply these results to the case of sufficiently small junc- 
tions, where we can assume Fj(A<p) to be a single-valued 
function of A<fi. Making a Fourier expansion of Fj(A<j)) 
and dropping terms based on the TR symmetry we get 

oo 

Fj{A(f) = -Ef - fiW cos(nA0), (2) 
n=l 

with a similar sine expansion for J s (A<f>). However, there 
exist "chiral" states, which are not TR invariant. An 
example is the A phase of superfluid 3 He, to which we 
return shortly. In such cases the TR symmetry no longer 
implies Fj(—A<f>) = Fj(A<fr) in general, but in special 
cases more general symmetries of the form Fj((3— A<f>) — 
Fj(f3 + A</>) may still be valid, JTT| see Section 5. 

An alternative expansion is to develop Fj in powers 
of the order parameters A L,R (similar to the Ginzburg- 
Landau expansion). Comparing this expansion with the 
Fourier expansion (0), we see that the coefficient Ej 
is proportional to A [1 + (3(A 2 )], where A denotes the 



amplitude of the order parameter. For small junctions 
the effective expansion parameter is the amplitude of the 
order parameter divided by temperature, A/fcsT c . Since 
this is small at temperatures close to T c , the lowest or- 
der terms dominate in the series (|J). Usually the first 
order term Fj = —E^p cos(A</>) is most important at 
all T, except in special cases where its amplitude is sup- 
pressed due to symmetry reasons. These symmetries for 
3 He or unconventional superconductors have been stud- 
ied in several papers, grecs 
of freedom of the order parameters but A<f> (and those 
related to the junction itself) are now embedded in the 
coefficients Ej, and these may be used as tuning pa- 
rameters. For example, in some cases[Q [l5| the sign of 

the (normally positive) E^p can be reversed, so that one 
obtains a "7r junction" , where the only minimum of a si- 
nusoidal EPR is at A(j> = ir rather than at A(f> — 0. In 
3 Hc-B such a trick may be done by controlling the order 
parameter textures with magnetic fields. H 

For vanishing E^ the higher order terms may still give 

(2) 

a finite critical current. In particular, a finite Ej term 
gives a 7r-periodic EPR due to the cos(2A^>) term. If the 
suppression of Ey is only partial, then some interesting 
mixtures of the 2n and 7r (and shorter) periodic com- 
ponents may be observed (see, for example, Ref. or 
Ref. ^|and references therein). In particular, it is exactly 
these higher-order terms of Eq. (||) which result in the 
"it states" of the p-wave (pinhole) junctions discussed in 
Refs. H ||, ^, 0. These tt states become more pronounced 
for T <C T c , since the higher-order Ey's are larger there. 
They are also the reason for the slanting of the CPR of 
an s-wave pinhole contact at low temperatures. |]l7f How- 
ever, in this case there is no way suppress E^ in order 
to single out the smaller tt periodic components. 

For large junctions we may expect two kinds of compli- 
cations to the description based on Eq. (|J). First, there 
is the effect of "kinetic inductance": as a result of cur- 
rent conservation, there must be a finite phase gradient 
carrying the current J s into, out of, and within the weak 
link.]l8| Thus the phase difference A<f> itself depends on 
J s , which makes Fj(A<f) and J s (A0) hysteretic (mul- 
tivalued) in general. As a result, the expansion in the 
cosines, Eq. (0), is no longer valid as such. The hysteresis 
can be modelled by introduction of the "slanting param- 
eter" and a self-consistent set of equations as in Refs. 
|l8| and |f9| Second, in the presence of a multicompo- 
nent order parameter the transmission properties of the 
weak link (reflected by the signs and absolute values of 

the coefficients E^p above) may depend on the order pa- 
rameter configuration inside the aperture. This configu- 
ration may change as a function of A<f>, which changes the 
CPR's.|| |9|, |2(| Thus the order parameter must be solved 
self-consistcntly in and around the hole, with boundary 
conditions specifying A L:R only somewhere far from the 
junction. The onset of these large- aperture effects in 
weak links of superfluid 3 He is the main subject of the 



3 



this paper. 

These changes can also be seen in the expansion of Fj 
in the order parameters. Instead of A/kgTc, the effec- 
tive expansion parameter in a junction of linear dimen- 
sion D turns out to be AD/Hvf, the gap divided by the 
(ballistic) Thoulcss energy. Using standard relations the 
expansion parameter can also be expressed as D/£gl(T), 
where £gl(T) is the Ginzburg-Landau temperature de- 
pendent coherence length. We see that the expansion 
breaks down when D / '£,gl(T) ~ 1. This is in agreement 
with our results below, which show multivalued EPR for 

III. SUPERFLUID 3 HE 

Below we shall only consider the triplet p-wave case 
where the gap vector A can be written as A„(x, k) = 

A M i(x)fci with a proper choice of the spin and orbital 
basis functions. pl|, E2[ In practice this means super- 
fluid 3 He, but a similar analysis may be valid in pos- 
sible triplet superconductors. In 3 He, in the absence of 
magnetic fields, there are two stable bulk phases, the A 
and the B phase. These are known to correspond to the 
ABM state A = A A d(xh + in) exp(i</>) and the BW state 
A = Asi?exp(i0), respectively. Here Aa,b ar e the bulk 
gaps of the A and B phases. In the B phase i?(o>, 9) is a 
rotation matrix, which may be parametrized by the rota- 
tion axis £: and rotation angle 9. In the A phase m _L n, 
and one usually defines a third unit vector 1 = rh x n 
which gives the direction of relative orbital angular mo- 
mentum of all Cooper pairs. Note that reversing the sign 
of d is equivalent to a phase shift by n, and that any phase 
shift by <fi is equivalent to a rotation of rh, n around 1 by 
angle —(f). The separation of the phase factor is therefore 
not unique, and it is further complicated by the existence 
of textures in the 1 field. However, in what follows, we 
can assume 1 to be constant most of the time, and if 
the same definitions are used consistently, no problems 
should arise. In the B phase the gap |A(k)| is isotropic, 
whereas in the A phase it has point nodes in the direction 
of the vector 1. To maximize the condensation energy, 1 
is therefore always rigidly oriented perpendicular to solid 
surfaces. The presence of a specific direction 1 for the 
Cooper pair orbital angular momentum means that the 
A phase is not time- reversal invariant. Consequently, in 
the context of weak links, some of the usually "obvious" 
symmetries must be reconsidered when A phase is in- 
volved. These symmetry properties are discussed below, 
in Section 5. 

The ABM and BW states are well-defined only in the 
hydrodynamic regime, i.e., on large length scales. Close 
to surfaces, for example, the order parameter will be 
modified on the coherence length scale £o due to scatter- 
ing of quasiparticles. A weak link involves length scales 
on the order of £o and a locally suppressed order parame- 
ter by definition, and therefore a more general treatment 



is required. In what follows, we let all the components 
of the order parameter vary freely close to the weak link, 
and the ABM or BW forms are fixed only on boundaries 
in the bulk liquid. 

A. Symmetries and Conservation Laws 

The free energy of a p-wave spin triplet Fermi super- 
fluid can be expressed as a functional of the order param- 
eter field A(x). In a region fi bounded by <9f2, we assume 
the free-energy expression to be of the form 

F a [A] = [ d 3 xf(A,VA). (3) 

Neglecting in this functional any terms that couple the 
spin and orbital degrees of freedom or introduce preferred 
spin directions (nuclear dipole-dipole and dipole-field in- 
teractions), it must remain invariant under global gauge 
transformations [E/(l)] and global spin rotations [50(3) s ] 
of A. These are parametrized by the independent real 
parameters <p an< i 9 a , a = x,y, z, respectively, which are 
independent of x. Infinitcsimally the transformations are 
written as 

A^ A'^ = e^A^ » A^i + iA^M (4) 
A^ > A^ — R^ v {59a)A U i ~ A^i — e a ^A U iS9 a , (5) 

where R(9 a ) is a (right-handed) rotation matrix. In this 
context Eq. (|^) should not be considered as a passive 
rotation of the spin-coordinate system, but rather as an 
active transformation of the physical state. Suppose that 
originally the order parameter corresponds to a station- 
ary point of Fq, i.e., it is a solution of the Euler-Lagrange 
equations 5Fq/SA — SFq/5A* = with fixed boundary 
conditions <L4|an = 0, for example. Then for variations 
8 A which leave / invariant, one finds a conservation law 
V-jjv = for some generalized current j^r. (This is a spe- 
cial case of Noether's theorem. P3[) In particular, for the 
variations of Eqs. @) and (j|) the corresponding currents 
are the "mass supercurrent" and the three independent 
components of "spin supercurrent" 

is = —{ +lA ^WA- +C - C j (6) 
df 

ja pm = +c.c, a = x,y,z. (7) 

Orbital rotations do not in general keep the energy in- 
variant and therefore no conserved "orbital supercurrent" 
exists. The physical interpretation of Eqs. (((]) and (0) as 
"mass" and "spin" currents (and thus their normaliza- 
tion constants) cannot be seen from this phenomenolog- 
ical approach, but they can be verified from microscopic 
theory. Note that the same current expressions can be 
obtained by inserting into Eq. (^|) the "gauge invariant 
derivative" prescriptions S^di — > S fJ _^d l +5^ v \.a i — e tiva b ai 
and expanding to linear order in a and b a i. 
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The conservation laws V • j s = and V ■ j^ pin = play 
an important role in the weak link problem since the cur- 
rent density distributions are closely related to the per- 
turbations of the order parameter at the junction. Due 
to the conservation of currents, the perturbations have, 
in principle, infinite range. In reality the symmetry- 
breaking dipole interaction makes spin currents decay on 
the dipole length scale £d. This is much larger than the 
scale considered in this paper and thus we can assume 
the conservation of both currents to be exact. 



B. Ginzburg-Landau Theory 

In numerical calculations we use the Ginzburg-Landau 
(GL) expansion of Fn, which is valid at temperatures 
close to the transition temperature T c , where the am- 
plitude of A is small. The GL theory and its ap- 
plications have been thoroughly discussed in various 
articles. [0 |l] [gj, [2^] Taking into account the bulk and 
gradient terms only, the GL energy density is 

f(A,VA) = 

— aTr(AA T *) + p 1 \Ti{AA T )\ 2 
+ /3 2 [Tr(yL4 T *)] 2 + f3 3 Tt(AA t A*A t *) 
+ 4 Tr(AA T *AA T *) + (3 5 Tr(AA T *A*A T ) 



(8) 



+ K&AlidjArt + K 2 diAl,diA 



+ KsdiA*. ,8^. 



This includes all the linearly independent terms up to 
fourth order (second order in gradients) which are invari- 
ant under global rotations of spin and orbital coordinates 
[SO(3) s '°], global gauge transformations [E/(l)], and un- 
der time reversal, i.e., complex conjugation. These in- 
clude the transformations of Eqs. (0) and (||). Other 
terms resulting from possible magnetic dipole-dipole 
and dipole-field interactions etc. could also be included, 
and they would not introduce major complications to 
the numerical calculation. However, these interactions 
only become important on length scales £d ~ 10 /xm 
and £,h which we assume to be much larger than the 
temperature-dependent coherence length £,gl{T) as de- 
fined shortly. This restricts the validity of our results 
from above on the field and temperature scales, since 
i H ~ H- 1 and S,gl{T) - (1 - TjT c )- x l 2 . In practice it 
may not be possible do accurate measurements so close to 
the critical temperature that the latter restriction would 
be a problem. Using the GL energy of Eq. (^) in Eqs. fl) 
and (u ), the formulas for the currents given in Refs. |24| 
and |25, for example, are exactly reproduced. 

The GL free-energy expansion was introduced above 
phenomenologically, with several parameters: a,f3i,Ki 
and 7. Values for these can be calculated from quasiclas- 
sical theory.^, ||(| All temperature-dependence of GL 
theory is in the coefficient a(T) = N(0)(1 - T/T c )/3, 
where N(0) = m*kF /2ir 2 h 2 is the normal-state density 
of states on the Fermi surface for one spin species. The 



gradient-energy parameters are 7 = 3, and K\/(^ — 
2) = K 2 = K 3 = K=(7((3)/2A0ir 2 )N{0)(hv F /k B T c ) 2 . 
These are all weak-coupling (WC) results, but the strong- 
coupling (SC) corrections will not be used, since they 
are not very accurately known, and are probably small. 
Also, K and a appear only in the natural length scale 
of GL theory, which we use as our unit of length. 
This is the temperature-dependent coherence length 
£ gl (T) = ^Kfc = ?(0)(1 - T/T c )-V*, where f(0) = 
-\/21£(3) /2A0TT 2 (hvp /k B T c ) is one way to define the zero- 
temperature coherence length £o- The WC values for 
the (3 parameters are 0f G /p B CS = (-1,2,2,2,-2), for 
i = 1, . . . , 5, where f3 BCS = (7C(3)/240tt 2 ) N(0)/(k B T c ) 2 . 
Pressure-dependent SC corrections to these have been 
calculated in Ref. [2?]. The main effect of the SC correc- 
tions is to change the difference A./ab = f%~ Fb °f ^ ne 
A and B phase condensation energies, f% = 2aA 2 1 /2 and 
f B = 3aA B /2, so that the A phase can become stable 
at pressures above the polycritical one, po- The value 
of po depends sensitively on the [3 parameters, and our 
fit gives it at roughly p = 28.7 bar, whereas experi- 
mentally po w 21 bar. Pressures close to po are impor- 
tant for studying weak links between the A and the B 
phases, since the phase boundary must remain pinned in 
the aperture. For a slit of width W the condition for the 
stability of the boundary is |A/ab| < 2(tab/W, where 
gab is the surface tension of the A-B interface. p8| 

The boundary conditions for the order parameter on 
solid surfaces generally satisfy A^Si = 0, where s is the 
surface normal. []l2| Everywhere in this paper we use a 
more strict boundary condition A 
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0. This corre- 
sponds to a microscopically rough surface, which scat- 
ters quasiparticles in a completely diffusive way. Most 
realistic surfaces are suspected to be of this type and, in 
addition, this is the simplest one to implement numeri- 
cally. 



IV. THE WEAK LINK PROBLEM 

We now apply the above considerations to describe a 
weak link of the form shown in Fig. |l|. It is a 2D approx- 
imation of the slit-type weak links used in experiments 
like those described in Refs. |^ and [l9|. This "archetype" 
is here considered between two infinite volumes, but we 
could also insert it between two flow channels or other 
restricted 2D geometries. [p9| The origin of coordinates 
is placed in the middle of the junction, with the x axis 
running through the aperture. Throughout, we call the 
sides with negative and positive x coordinate the left (L) 
and right (R) sides, respectively. 

The region Q in Eq. (||) is now the one inside the 
outer circular arcs of radii i?oo. In 3D the i?oo cutoff 
could be taken to infinity. In 2D an inconvenient finite 
value for i?oo is required due to the logarithmic ln(i? oc /r) 
(and n ot 1 /r) as ympt otic behavior of the phase fields [see 
Eqs. (A5) and (A13) in Appendix A]. This cutoff may 
be thought to describe some effective length scale, deter- 
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FIG. 1: A 2D representation of the slit-like weak link in xy 
plane. The circular arcs are the cutoffs used in the numerical 
calculation. is a necessary artifact of the 2D approxima- 
tion, and it is on the order of the z directional length L of the 
assumed slit-type weak link. Q In a 3D calculation one could 
choose Roo = oo without problems. R c is an artificial cutoff 
dividing the computational effort into numerical (inner) and 
analytical asymptotic (outer) regions. In the A phase, where 
the 1 vector will be perpendicular to the wall, the circles will 
be compressed to ellipsoids with the scaling of the perpendic- 
ular coordinate x by v2- 



mined by the dimensions of the container, or a distance 
after which the decay turns three-dimensional (1/r), for 
example. The z directional length of the assumed slit 
gives a good approximation for the latter. |J In the A 
phase a more natural choice for the form of the cutoff 
arcs is actually ellipsoidal, with the perpendicular coor- 
dinate scaled down by v(7 + l)/2 » V2. This follows 
from the smaller superfluid density in the direction of the 
orbital angular momentum vector 1 — see Appendix A. 



A. Coupling Energy and Soft Variables 

The left and right boundary conditions are in general 
functions Aj^^sl), A^sr) of sl,r which parametrize po- 
sitions on the T L,R . The Josephson coupling energy is a 
functional Fj[A L , A R ]= min^FoLA] — Fq. Here we ex- 
plicitly subtract a term Fq, which is the energy in f2 for 
the same A L ' R , but with the aperture closed. This is 
done to remove the bulk and surface contributions, but 
in 2D Fj still depends on R^ due to the currents, as ex- 
plained in Appendix A. Here we have assumed that fixed 
boundary conditions are suitable, and that they define 
the minimum uniquely. Actually, there can be several 
local minima which are separated from each other by en- 
ergy barriers. Therefore the current-phase relation may 
be multivalued, each CPR branch corresponding to a dif- 
ferent type of minimum. Jumps ("phase slips") between 
these branches are hysteretic. 

As special cases we shall consider the A and B phases 
of 3 He. There are then four main cases which can re- 



alized for a hole in a nonmagnetic planar wall. Using 
pairs of letters to denote the L and R phases, respec- 
tively, these are BB, AA with parallel (TT), AA with 
antiparallel (fl) l's, and the different configurations of 
AB. The asymptotic forms in all of these cases can be 
treated very similarly by writing A^i = R^A^ (a;)e 1( ^, 
where the broken-symmetry variables, or "soft degrees of 
freedom" R((2>,0) and <fi are assumed to be constants on 

The function A^ (x) is the order parameter calculated 
for a planar diffusive wall in the absence of a weak link, 
and thus includes the suppression at the wall. A one- 
dimensional minimization is required to determine it (see 
Appendix B). In the B phase A^ = diag(A_L, A||, A||), 

oo) — ArSm in bulk. In the A phase 



with A[ ){x 



,(°) 



m|°h, where rrS ^ 



AO) 



A °i = a{x)d ux [ml 
±Si Z and a(x — ► oo) = A^. Here the A and B phase bulk 
gaps A a and A^ are defined in Appendix B. The positive 
(negative) sign of nf^ chooses the vector I = m x n to 
point in the positive (negative) x direction, which are two 
degenerate configurations. A more familiar form for the 
spin part of the A phase order parameter is obtained by 
defining the d vector with = R^ V 8 UX = R^x- 

Assuming A^ to be fixed, it suffices to write Fj only 
in terms of the bulk order parameters, or more precisely, 
the soft degrees of freedom <fr and R. If we require 
Fj[A L ,A R ] to be invariant under global gauge trans- 
formations and spin rotations, it should be unchanged 
when A L ' R are both multiplied by R L e~'^ . Thus we 
have Fj = F J [A^ L (x),R L ' 1 R R A^ R (x)e iA< l'], where 
A<p = (j> R - 4> L . Now, since A^ L and A^ R are fixed 
by assumption, we see that Fj can be parametrized sim- 
ply as 



Fj=F J (Acf > ,R%Ri 



h3 =x,y,z 



(9) 



- uL 



R L R R 



where an orbital-space rotation matrix ib«j 
i,j = x,y,z remains as an argument. In Ref. H this was 
derived only for the BB case, but we now see that it is 
valid more generally. For the A phase only the three com- 
ponents = R^x out of the nine R^i are relevant. This 
implies that for the AA case Eq. (0) can be simplified as 



Fj = Fj{A<j>,d L -d R ) 



and for AB 



Fj = Fj(A^O), 



(10) 
(11) 



where we have defined an orbital-space vector Oi = ip x i = 

Often it is also instructive to consider the following 
"tunneling" form of the coupling energy: jl^l 

FT = Re[aA^ x A R x + bA%A% + cA^ z A R z }. (12) 

This is the lowest order term (oc A 2 ) in the order- 
parameter expansion discussed in Sec. 2. Thus it is also 
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the leading term in the Fourier expansion 



F 



(i) 



-E { p cos 



^tun 



(9(A 4 ). There are a few more 



assumptions coming into Eq. ([l2j). The suppression of 
the order parameter near walls causes an ambiguity what 
location r should be used for A(r) in Eq. (|l^). Here we 
have used A L ' R in the bulk, and defined A = A/A^.b de- 
pending on the phase. Also, we are assuming orthorhom- 
bic symmetry [(2/m)(2/m)(2/m)] as for the slit above. 
The coefficients a, b and c are different for the different 
phase combinations and temperatures, and c = b if the 
symmetry of the junction is cylindrical [(oo/m)(2/m)] 
instead. Q 

B. Josephson Currents 

Assume now the situation of Fig. either 2D or 3D. 
There can be no current flow through solid walls, and 
thus j s • s = j^ pm • s = on surfaces, where s is the local 
surface normal. Therefore, the conserved DC mass and 
spin currents flowing through the junction are given by 
(dS • f >{)) 



Js = 



rspin 

J rv 



dS 



•spin 



ds-j, = - / dS-j s (13) 

dS-j s / n - (14) 

These inherit their symmetry properties from the Joseph- 
son coupling Fj, and the expressions revealing this de- 
pendence can be derived as follows. 

Consider infinitesimal variations of i*h[A] around the 
stationary point, which are of the form shown in Eqs. (Q) 
and (^) but now with local variational parameters <5</>(x) 
and 59 a (x). We choose them so that 54\t l — 59 a \ T L = 
and 5cj)\ r R = Scj) R , 59 a \ r n = S9 R . Due to the stationarity, 
only the boundary term contributes to linear order, and 
using Eqs. @ and (0) we have 

(15) 



SF n = (h/2m 3 )JM R ~ JT n S8 R . 



On the other hand, we can similarly transform the 
boundary values A</> and ipij, so that SAtfi = 5tfi R 
and Sipij = -e af3l R^R R 59 R . By expanding Fj(A(j) + 
5(j) R ,ipij + Sipij) — Fj(A(f),ipij) + SFj to first order in 
5<j) R , 56~ 



J, 



and equating SFj with Eq. (|15|), we find 
2m 3 dFj 



H dA<t> 



dF 



J 



rspin r>L r>R 



(16) 

a = x,y, z. (17) 



For a BB junction these expressions were derived in Ref. 

but again we see that they are more general. Using 
the simplified forms of the coupling energy for AA [Eq. 
©] and AB [Eq. (0)] cases, Eq. @ reduces for an AA 
junction to 



jspin = x d flj c 



dFj 



d(d L ■ d R ) 



x,y,z, (18) 



and for an AB junction and to 



Tspin 



£a/3-tdpR~[i 



dFj 



a 



x,y,z. 



(19) 



Note that the expansions of Fj leading to Eqs. ( p^ ) 
and ([[7]) work at most locally, i.e., on each branch of 
the solution-space separately. In some simple cases one 
might also use them in order to check the consistency 
of one's numerics. To remove the ambiguity related to 
the A-phase phase factor, we shall always choose 1 = ±x 
with n(°) = ±z, as explained above. This is a natu- 
ral choice since n-reversed configurations are related by 
time reversal. When the phase factor is included, TR 
should be understood as a reversal of both <j> and n^ ^. 
On the other hand, a reversal of n(°) alone is equivalent 
to reversal of 1 at constant 6. 



C. Controlling the Boundary Conditions 

Above we assumed the order parameter at the Roo cut- 
off to be fixed to the form A(x) = RA^(x)e i,p . Here A®) 
determines the bulk phase, including its modification at 
walls. In case of the A phase, A^ must always be such 
that 1 || ±s, where s = ±x is the wall normal, to minimize 
loss of condensation energy. Besides the phase angle, this 
leaves only the rotation matrix R to be determined by 
some hydrodynamic interactions. In the A phase this 
reduces to fixing d with a combination of the dipolc- 
dipole oc — (d • l) 2 and the dipole- field oc (d • H) 2 bulk 
terms. Unfortunately, for H || 1 the configuration remains 
undetermined. Furthermore, it may be difficult to pro- 
duce the most interesting situations where d L x d 7^ 
by any physical means. In principle, one way is to use 
magnetic fields of different directions on the two sides. 
(Some more exotic ways, such as an A-B interface, could 
be imagined. ||) In the B phase there is more freedom 
in controlling R. To start with, we always assume the 
rotation angle 9 to be in the minimum 9q 0.587T of 
the dipole-dipole oc (cos 9 + 1/4) 2 interaction far from 
the junction. |3[ The remaining uj vector is coupled to the 
wall (with normal s) by a surface-dipole term oc — (u>-s) 2 , 
and to magnetic field via a bulk term oc — (lj ■ H) 2 and 
a surface term oc — (H • i?s) 2 . As discussed in Ref. ||, 
it is the surface-dipole and the surface-field terms which 
determine Cj close to surfaces at low and high fields, re- 
spectively. For the magnetic surface configurations we 
use the definitions A-D of Ref. ^, but refer to them with 
lowercase letters a-d. 

In addition to the rotation matrix, also the phase dif- 
ference A(j) must be specified. In a channel geometry 
this would be replaced by specifying the phase gradient, 
i.e., the superfluid velocity^ We assume that the in- 
verse Josephson frequency oj^ = H/2A/J, is much larger 
than all order parameter relaxation times (except that of 
Ac/)). In practice this should be well satisfied, since ujj 
tends to be no higher than in the audio regime. Q This 
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makes the equilibrium concepts of energy-phase Fj(A(f>) 
and current-phase J s (A<p) relations sensible, and there- 
fore warrants the present calculation. However, the situ- 
ation in the A phase may again be more complicated due 
to the orbital viscosity phenomenon, which slows down 
the dynamics of the I vector.^) Therefore we cannot ex- 
pect the calculation to properly describe the dynamics of 
phase slips, i.e., jumps between branches of J s (A(j>). 



spin x current 





D. About the Numerical Implementation 

The cutoff Roo is arbitrary, but from the CPR's for 
one choice, we may determine the CPR's for any other 
one.|)| For example, we may do the numerical calculation 
for r < R c where R c i s small (see Fig. |l|), and then use 
Eq. (A7) or Eq. ( Al5| ) to correct the phase differences 
for some r — Roo ^> R c , so that the new CPR becomes a 
"slanted" version of the original. Whenever spin currents 
are present, the spin-rotation boundary conditions should 
also be modified using Eq. (|A|) or Eq. ( |Al~6| ). However, 
in the B phase this will lead to rotation matrices which 
are not of the form R(uj,9q) that we wish to have for 
large R^, (although still Roo < £d). 

Thus the calculation should in general be done for each 
(large) cutoff radius Roo and each set of spin-orbit bound- 
ary conditions separately. The numerical minimization 
may still be done inside R c <C Roo, and for r > R c the 
asymptotic solutions of Appendix A are used. In our cal- 
culations, the fitting of the solutions at R c is not done 
by comparing derivatives at any location on the bound- 
ary, but by the more physical requirement of conservation 
of total mass and spin currents over R c . This leads to 
a somewhat cumbersome self-consistent iteration proce- 
dure, which is described in Appendix B with the rest of 
the numerics. Fortunately, these practical issues should 
not affect any of the above or the following analysis, ex- 
cept through the value of Roo, and thus a discussion of 
the solutions in the asymptotic region is postponed to 
Appendix A. If these "asymptotic corrections" are not 
used, then Roa — R r - 



V. CPR'S FOR AA, BB AND AB JUNCTIONS 

In this section, we present a more careful analysis of 
the BB, AA, and AB junctions. In what follows, we 
always present results for BB at vapor pressure p = 
bar, for A A at melting pressure p — 34.4 bar, and for 
AB at roughly the coexistence pressure of p = 28.7 bar 
(Af AB = 0). 

In numerical calculations of the CPR's, the asymp- 
totic corrections were usually taken into account, and 
we used the outer cutoff Roo/S,gl ~ 30. In this case 
the value for the inner cutoff was roughly R c /£,gl = 10, 
although in principle the results should be quite inde- 
pendent of it. As expected, the logarithmic decay of the 
phase corrections is slow, and as we varied R c /£,gl and 
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FIG. 2: Mass and spin current density distributions for BB 
with W/^gl = 8, D/^gl = 4, -ci L = u R = x, and = 
0.3-7T. The surface spin currents of B phase are clearly visible. 
Close to edges, small spin whirlpools are often formed. 



Rool^GL m the range 10 . . . 60, only small differences in 
the CPR's could be seen. Most of the time a lattice spac- 
ing of Ax/^gl — Ay/^GL = 0.5 was used. Refinement 
did not lead to qualitative differences in the form of the 
CPR's, and only to differences of at most a few percent in 
the critical currents. The minimization was usually car- 
ried out until the error (as measured by the norm of the 
gradient in the R c region, see Appendix B) was smaller 
than 10 -5 . For given boundary conditions, an accuracy 
better by more than 10 orders was possible, but com- 
pletely unnecessary, considering the uncertainties related 
to the handling of the bulk cutoffs. Current conserva- 
tion at every lattice point was practically exact, except 
for points on the R c cutoff, where only an accuracy of 
|V • j s | ~ 10~ 2 was achieved locally. However, the total 
current over the cutoff was required to be conserved down 
to 10~ 4 , if the asymptotic corrections were used. These 
values refer to the p s and p± units used in Appendix A 
and B, which are used also in all the figures below. 



A. BB Junctions 

The case of BB junctions was already considered in 
Ref. [jj, and we skip most of the analysis. Note that for 
BB time-reversal symmetry (complex conjugation of A) 
implies the simple relation Fj(-Acj)) = Fj(A(j>) in ad- 
dition to 2ir periodicity, and thus J s (—Acj)) — — J s (Acj)), 
J^ in (-A0) = JJ in (A0). Inserting Eq. Q with A = 
i?exp(i</>) into Eq. (O) gives the spin current as J„ pm = 



w [aR 



L R R 

fix ^VX 



] cos Acj). Notice 
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FIG. 3: BB current-phase relations for W/^gl = 2, D/£gl = 
2, uj r — x, and uj l in xy plane with azimuthal angle r/L 
with respect to x axis. The spin current components are as 
follows: x (open circle), y (plus), z (closed circle). Note that 
there are two scales on the bottom axis: a different CPR with 
A(/> in range . . . 2-7T is shown for each of the values t/l/tt — 
0,0.2,0.4, ...,1.0. 
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FIG. 4: BB current-phase relations for W/^gl = 4, D/£gl 
2 - otherwise same as in Fig. |^. 



that even if Fj 1 ^ is identically zero, and thus the EPR is 
7T periodic, the spin current may remain 2ir periodic. 

Despite the apparent simplicity, the BB junction has a 
rich spin current structure. This is because, in addition 
to the Josephson spin currents induced by bulk boundary 
conditions, there exist spontaneous spin currents close to 
solid surfaces. |n] An example of these is shown in Fig. ||, 
where small spin current loops are seen to be stabilized 
in different parts of the weak link. However, these are 
not "topological" vortices with a quantized spin velocity 
circulation driving them. The cores are usually associ- 
ated with a small A-phase-like "orbital magnetization" 
oc eijk lm(A* j A flk )/2 along the z axis. 

Figures || and |J show the CPR's for apertures of two 
different widths W, but the same "depth" D, i.e., thick- 
ness of wall, using the same series of boundary conditions: 
CJ R = x, and ui is in the xy plane with azimuthal an- 



Im A 



FIG 



5: Schematic x dependence of Im A for k = x at a> = 
u>" (9 — 0) and A</> = n on the branch (a) and n branch 
(b). For both cases Re A = 0. With configuration (b) a lower 
energy may be achieved, since it avoids a singularity where 
the gap vanishes. 



gle r]L from the positive x direction. These figures are 
equivalent to Fig. 11 of Ref. |5| and, while corresponding 
to rather imaginary experimental conditions, they show 
quite clearly the essential differences between the CPR's 
of narrow and wide apertures. In the narrow case (Fig. 
||) the behavior follows well the perturbative picture of 
Eq. (0): the energy-phase and current-phase (mass or 
spin) relations are non-hysteretic and close to 27r periodic 
(co)sine functions. Only in cases where Fj un oc aR^ x R R 

i ldL D-R 



cR^ z R R z vanishes (as it obviously does for 
r]L ~ 0.4-7T and t)l ~ 0.77r) a small 7r periodic contribu- 
tion due to the remaining higher-order terms is visible, 
but with very small critical current. These higher har- 
monics become stronger and produce more pronounced 7r 
states at low temperatures.^] For r] L w 0.67T the CPR is 
sinusoidal but u tt shifted" , whereas for rj L w ir the ir shift 
again vanishes. This is different from the pinhole case of 
Ref. |j| where a ir shift is obtained for uj l — uj r = ±x. 
This is apparently due to the difference in the symme- 
try of the aperture: unlike for the pinhole, here we have 
c^b. 

On the other hand, in the wider case (Fig. |j) the CPR's 
are already clearly hysteretic. The leftmost panels have 
CJ L = lj r = x, which is the situation studied in Ref. || 
for the circular 3D aperture. Here the behavior is simi- 
lar: in addition to the usual "0 branch" there exists also 
a "7r branch" around A(j> = ix. Figure ^| depicts the corre- 
sponding behavior of the order parameter inside the aper- 
ture at exactly A(f> — ir (compare with Fig. 3 of Ref. |^). 
Here the constant rotation matrix has been dropped (i.e., 
9 = 0) so that A is diagonal at the infinities. On the "0 
branch" the behavior of A(x, k) is singular, whereas on 
the 'V branch" A avoids a singularity by "escaping into 
the third dimension" . The rotation of A around y with 
growing x means that there is a spontaneous Josephson 
current of y directional spin flowing through the aper- 
ture. Adding the spin rotation by 9 = 0.587T around 
x gives the y and z spin current components visible in 
Fig. |J. As pointed out in Ref. 0, the order parameter on 
the 7T branch has a similar structure as in a double-core 
vortex, |25| where the virtual vortex cores are inside the 
wall, one on each side of the channel. 

For parallel uj vectors there can never be a jump from 
the to the 7r branch, and so the conclusion in Ref. || 
was that the branch may not be experimentally achiev- 
able. But this is not so. Once the rotational symmetry 
of the bulk boundary conditions around the surface nor- 
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mal is broken by perturbing Cb even slightly, the hys- 
teresis of the branch is reduced, and the tt branch be- 
comes the branch of minimum energy for A0 w tt. Since 
J' s (A<p = tt) > 0, this "tt state" is even a (local or global) 
minimum of EPR and thus stable against small pertur- 
bations of A</>. Notice how in Fig. || the order parameter 
transitions between and tt minima continue to be as- 
sociated with changes in the spin currents, although a 
simple interpretation as with Fig. |5| is no longer possible 
for large tilting angles of u) L . But for small angles a jump 
onto, or away from the tt branch can be interpreted as 
a phase slip by a half-quantum vortex.^ Note also that 
A<p = tends to remain at least a local minimum of 
the EPR, i.e., J' s (A(f> = 0) > 0, so that a sinusoidal 
but apparently tt shifted EPR ("7r-j unction") tends to 
be avoided. In the narrow-aperture limit W/£gl ^ 3 
the situation J' s (0) < appears to be more easily ob- 
tained (compare Figs. || and ||). The last panels of Fig. |] 
correspond to — Cj l = uj r = x, for which a tt branch 
with J' s {Atf> = tt) < was found already in Ref. ||. 
Our simulations confirm this old result, and show that 
the J' s (A(f) = tt) < branch of Ref. || is related to the 
branches with J' s (Acf) = tt) > (stable tt state) by a 
continuous variation of the boundary conditions. 

The limit between the two types of behavior, Figs. || 
and |, is rather clear-cut, and it only depends on the 
width W. The transition occurs roughly at the width 
W/£gl = 3, which coincides with the analytic criti- 
cal value tt for destruction of superfluidity in an infinite 
slab.|32|, [3^] However, in the weak link case the order pa- 
rameter amplitudes and critical currents always remain 
finite, although small. In wide apertures the CPR's with 
tt states can have relatively large critical currents — com- 
parable to those of branches in general. In the narrow 
apertures, however, there tends to be almost an order-of- 
magnitude difference between E^ 1 ' and E^> in the GL 
regime. 

Figure || shows the CPR's for a more physical control 
parameter, namely the angle Oh of an applied magnetic 
field. Here we must assume that the field is strong enough 
for the texture to be determined by the magnetic surface 
interaction — (H • i?s) 2 , but small enough for ^> £gl, 
since we have neglected the magnetic terms from our GL 
free energy. We choose to show the "bd" configuration 
defined in Ref. ^ with H in zy plane, since this is an 
example of a "tt junction" with W/^gl > 3. In fact, a 
rather similar figure as Fig. |^ exists also for the "ad", 
"ac" and "be" configurations with field in the plane of 
the wall. For the "ab" or "cd" configurations no tt states 
nor tt shifts were seen, although many magnetic field di- 
rections (both in-plane and off-plane) were checked. 

Based on these results we conclude that a rich variety 
of tt states and tt shifts can exist in single BB apertures 
in the GL regime, although it may be difficult to realize 
the required Cj configurations experimentally. Approxi- 
mately the behavior of Fig. where the unreachable tt 
branch becomes visible, could be obtainable by slowly 
turning on H while keeping it in a suitably chosen direc- 
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FIG. 6: BB current-phase relations for W/^gl = 4, 
D/^gl = 2 in the magnetic "bd" configuration, field in 
the yz plane with polar angle 8h = 0.0, . . . , 0.5-7T. For 
8 H = thus u L = (-^/T/5,-^/3/5,-^/T/E) and uj r = 
(— \/l/5, +\/3/5, +\/l/5) (see Ref. ^, but note the permuted 
coordinate axes). 



tion. Without a doubt, a qualitatively similar division 
into small and large aperture behavior will be valid at low 
temperatures also. To verify this, a quasiclassical calcu- 
lation for the general-size weak link would be required. 
However, this appears to be too demanding considering 
that there is no reason to expect any essentially new phe- 
nomena. The critical currents and a phase diagram for 
BB are summarized in the next section. 



B. A AfX and AAfj Junctions 

In AA junctions the time reversal symmetry is no 
longer quite so simple as for BB, due to the fact that 
complex conjugation also reverses the direction of 1. In 
both the AA| { and A Aft configurations TR reduces only 
to Fj{\ L ,l R ,(3 - A(f>) = Fj(-l L ,-l R ,(3 + A(j>), where 
(3 = or tt. (Note again that reversing 1 is done by 
flipping the reference direction n' '.) However, revers- 
ing the directions of both l's simultaneously has no ef- 
fect in our orthorhombically symmetric junction, so that 
Fj(-l L ,-l R , A<t>) = Fj(l L ,l R ,A(t>), and therefore the 
TR symmetry reduces to that of the BB case. The 
current-phase relations therefore satisfy J s (f3 — A(j>) = 
- J s (/3 + A<j>) and J^ pin (/3 - Acj>) = J^ pin (/3 + A<f>) as for 
BB. Since ±1 L ' R correspond to the same branch(es), the 
symmetry also implies that the minimum of EPR is at 
either Aip = or tt, depending on the directions of the d 
vectors. In our case A<fi = is the minimum if d = 1 for 
both L and R. 

For AA|| the critical current will aways vanish if the 
aperture has full rotation symmetry around the wall nor- 
mal s = ±x. This is because an orbital rotation by 
9 around s adds a phase =F# to the order parameter if 
1 = ±s, and therefore A<fi = <fi R — <f> L is changed by ±28. 
Since this rotation is a symmetry operation for the aper- 
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FIG. 7: Different configurations for the 1 field in an aperture, 
(a-c) are for antiparallel and (d-f) are for parallel l's. 
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FIG. 8: The top panels show the 1 texture for W/^gl = 6, 
D/£,gl = 4, A<j> = 0.3tt for AA U (left) and the same for 
AA 11 (right). The corresponding mass currents are shown 
in the bottom panels. 



ture, it should not change Fj ([H]). Thus Fj should have 
the periodicity Fj(Acj)+29) = Fj{A<j>) for all 9. This can 
only be satisfied if the phase dependence of Fj vanishes 
altogether. In general, if the aperture has n-fold rotation 
symmetry (n > 1), the fj, EPR's are 2(27r/n) periodic. 
For example, for the orthorhombic slit considered in this 
paper n = 2 and only a rotation by 9 = Tt must be a sym- 
metry operation. This implies no additional restrictions 
on Fj, since it is already 2ir periodic. The case AA|| has 
nontrivial 2ir periodic Fj's irrespective of the geometry, 
except for the special case d L ■ d R = (see below). 

Figure ^| represents the different configurations for 1 in- 
side the aperture. The middle one for both f J, (b) and |T 
(e) is symmetrical with respect to the x axis. Case (b) is 
very singular and does not usually correspond to the min- 
imum energy in large apertures, but may be metastable 
at least for A<fi ss 0. Here a radial singularity is depicted, 
but a "hyperbolic" one is also possible. In some cases 1 
can also escape from the plane. The H configuration (e), 
on the other hand, is usually the minimum-energy con- 
figuration for W 3> D and A(f> w 0, as may be expected. 
The other configurations correspond to pairs of degen- 
erate cases, where 1 bends asymmetrically, as shown in 
Fig. H in more detail. Actually, the quantity shown is 
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FIG. 9: AAf-i current-phase relations for aspect ratio W : 



D = 2 : 1, 1 — d = —x, 1 = 
scaling parameter I = D/(4£gl)- 



x, as a function of the 
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FIG. 10: AATT CPR's for W/£gl = 6 and d L 



L,R 



as a function of D/^gl- For small D, the symmetric 1 texture 
is the minimum-energy configuration, whereas for large D it 
is one of the asymmetrical ones. In the range D/^gl — 2 ... 4 
the configuration depends on A<j>, which makes the CPR dis- 
continuous. Note that for D/^gl > 3 the critical current no 
longer changes, but the CPR becomes slowly more hysteretic. 
Compare with Fig. [l4| below. 



eijk Im(j4* ■j4 /1 fe)/2, where A = A/Aa, which reduces to 

li in the bulk. Figure ^ also shows the mass current den- 
sity fields associated with a phase difference A(j> = 0.3tt. 
Because the superfluid density is twice as large in direc- 
tions perpendicular to 1 than for those in parallel with it, 
the bending of the 1 fields results in asymmetrical current 
distributions close to the junction. These effects result 
in a more complicated numerical problem than for the 
isotropic B phase (see Appendix A). 

Figure || shows CPR's for H and an aperture of fixed 
aspect ratio, but different dimensions. No ir branches ex- 
ist as for BB, and the corresponding figure for ff would 
be very similar. It was noticed that, both in fj, and 
|f cases, when a jump to another branch of J s (A</>) oc- 
curred, the asymmetric texture shifted from one degen- 
erate configuration (a or d) in Fig. [7] to the other (c or f). 
At least in the || case this could be understood by re- 
membering that a phase slip of A<f> by 2tt is equivalent to 
a rotation of the orbital space around s by tt. However, 
these results are due to the minimization algorithm, and 
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do not necessarily correspond to physical dynamics. 

As pointed out above, for A A "ft at small D and 
W/^glJ^j 3 the minimum-energy configuration is that 
of Fig. 0e when Atfi = and no current flows. But when 
A^ is increased, a change to one of the asymmetric 1 
configurations (d or f) will occur, which causes a sudden 
reduction of the current and thus a discontinuity in CPR. 
This is illustrated in Fig. |To| — see also Figs. [14| and |l5| 
below. The larger D is, the smaller is the critical A<fi for 
this event. An accurate calculation of this critical value 
is difficult, however, since it depends on whether or not 
the asymptotic correction for 1 is used, and on how far 
the R c cutoff is taken. This introduces some amount of 
additional uncertainty in the forms of the CPR's. For 
AAfJ, one of the asymmetric configurations (a or c) is 
always the minimum one (except for small D and W , see 
Fig. |l5|) and thus in Fig. |] all the CPR's are smooth and 
well-behaved. 

For d L • = the CPR's are tt periodic in both AA 
configurations, as can be expected by noticing that Fj un 
[Eq. (p^[)] vanishes in this case. However, as discussed in 
Section 4, there appears to be no easy physical way in 
which a situation with d L x d R could be achieved. 
According to Eq. (|l8| ) this more-or-less rules out spin 
currents, and thus seems to render the AA CPR's rather 
un-interesting, apart from the effects of Fig. [l^ arising 
from the 1 texture. Theoretically this is not a problem, 
of course, and if we force the configuration d L • d^ = 0, 
then the following observations can be made. In the 
small-hole limit (CPR's continuous) the critical currents 
are considerably suppressed, and J s (A<f)) w J c sin(2A0). 
However, in the large- hole limit (CPR's hysteretic), the 
critical currents are not significantly affected. The CPR's 
look just like in Figs. ^| or [h], but with the tt branches 
being added in between each branch. As the d's are 
gradually changed from d L ■ d R = 0, either the or the 
tt branches begin to rise in energy, and the correspond- 
ing CPR branches get smaller (similar to metastable "tt 
states" in BB). Finally the branches become unreach- 
able and 2tt periodicity without any tt states is regained. 
Again, these observations hold for both AA"f J, and AA|| 
alike. The most notable difference between them is in 
the critical currents, which are summarized in the next 
section. 



C. AB Junctions 

The case of an AB junction is the most complicated 
one to analyze. In this case the requirement of time re- 
versal symmetry implies Fj(l, j3 — A(f>) = Fj(— 1, j3+ A<f>), 
where (3 = or tt as in AA. This no longer simplifies to 
the intuitive BB form as it did for our AA junctions. 
Rather the EPR's for 1 = ±x form completely separate, 
unsymmetrical branches, which are 2tt periodic in gen- 
eral, and which can both still separate into hysteretic 
sub-branches as usual. The minimum of the EPR is also 



not restricted to Acj> = or tt on either branch (see Fig. 
[Tl| below). Similarly, the currents only satisfy relations 
between the ±1 branches J s (l, 0- A(j>) = -J a (-1, /3+A0) 
and J^ pin (i,/3 - Acj)) = J^ in (-l,(3 + A<j>). The roles of 
A</> = and tt can again be interchanged by flipping the 
d of the A phase. 

To analyze the AB case further, we consider the tun- 
neling term, Eq. (p"2|). Assuming A phase on the left 
and B phase on the right, we again introduce the or- 
bital space vector Oi = d^R^ and find that Fj un = 
Re{{bO y (m y ± ih y ) +cO z (rh z ± in z )] exp(iA^)}. We see 
that if dfj, = zLR^iSi, i.e., O || s, then O _L m, n and 
Fj nn vanishes. This is true, for example, in the most 
symmetric (and, in the absence of magnetic fields, most 
probable) situation where both d and the B phase rota- 
tion axis u) are perpendicular to the wall, but not neces- 
sarily if either of them deviates from this configuration. 
Incidentally, — ±R^Si is also the minimum-energy 
configuration for a free, planar AB interface, when s is 
the interface normal. [Q 

More general conclusions can be based on the general 
form Fj (|ll| ) in the symmetric case 6 || s. Similarly as in 
the case of AA||, a rotation of the orbital space around 
s by angle 8 is now equivalent to a shift of the phase 
difference <j> by ±6>. (The factor is different here since the 
B phase side is not affected by the rotation.) In the case 
of a circularly symmetrical aperture we must then have 
Fj(Acf> + 8) = Fj(Acf)) for all 8 which, again, means that 
the phase dependence must vanish to all orders. For the 
slit junction this implies only Fj{A4>+tt) = Fj(Acj)), i.e., 
that the energy-phase relations are tt periodic. Generally, 
if the aperture has n-fold rotation symmetry, then the 
CPR's must be 2tt/ti periodic. However, these simple 
conclusions no longer hold when O | s, and in general 
the nontrivial 2tt periodic behavior is obtained in any 
aperture geometry. 

Figure [ll] shows examples of CPR's when d and lJ vec- 
tors are controlled by a field H strong enough so that d 
points in the direction corresponding to minimum pos- 
sible dipole energy — (d • l) 2 allowed by a strict condi- 
tion d ■ H = 0. Here Cj is chosen to be in the "a" 
configuration. H The failure of the phase- inversion sym- 
metry Fj(—A(f>) = Fj(A</)) is clearly visible; Figure [ll] is 
for 1 = x, and the branch for 1 = — x is obtained by us- 
ing the above symmetries for Fj . Figure [fiL on the other 
hand, shows the CPR's for a fixed aspect ratio 2 : 1 (as in 
Fig.|| for AA||) but different dimensions in the 6 || ±s 
case. The CPR's are tt periodic as they should, and even 
show the presence of strong 7r/2-periodic admixtures, or 
separate u tt/2 branches". 

The order parameter profile on the x axis which goes 
through the AB interface, and the form of the 1 field 
there are shown in Fig. (l3l Here we have A Jab > 
(p > Po) f° r illustration purposes, although the CPR's 
were calculated for Aab — 0. The boundary of the 1 
field, i.e., the A-B interface is thus seen to bulge into the 
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FIG. 11: AB current-phase relations for 1 L = x, W/£gl = 4, 
D/^gl = 2 and B phase in the magnetic "a" configuration,^] 
with field in the yz plane with polar angle 8h- Thus d = 1 
for all 8h, and, for example, the B p hase rot ation axe s for 
H = and 6 H = 0.5tt are w = (+ y/IJE, - s/S/E, 
and ui = (+y/T/5, +s/l/5, + ^/3/5), respectively. The EPR's 
and spin CPR's for l L — — x branch (not shown) could be 
obtained by mirroring with respect to A<f> = n, and the mass 
CPR's with the mirroring and a change of sign (see text). 
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FIG. 13: The 1 field terminating at the pinned A-B interface 
for W/£gl = 20, D/^gl = 4 at p = 34.4 bar. Also shown is a 
slice (at y — 0) of the order parameter through the interface. 
Here Atf) = and d M = R^iSi, where the constant R matrix 
has been dropped in both phases for clarity. The cut-off is at 
Rc/^gl = 30, i.e., at x/^gl = — SO/y/2 in the A phase. The 
asymptotic corrections make the order parameter transition 
at the cutoff smooth. 



symmetry with respect to the xz plane and 1 pointing 
out of the plane in the ±z direction at the interface were 
also often seen on some branches. 



VI. SUMMARY AND PHASE DIAGRAMS 



FIG. 12: AB current-phase relations for aspect ratio W : D — 
2 : 1, and \ L — d L = +x, lj r = x (or = ±R^,iSi in general) 
as a function of the scaling parameter I — D / [A^gl). 



B phase. Even at A/ab = the interface always settles 
into the B-phase end of the channel. p8[ As in the AA 
cases, 1 tends to bend parallel to y inside the hole. This 
is the preferred configuration also for the A-B interface 
itself. The order parameter components shown on the 
right are the same as in Ref. [3^ or Ref. [36[ except for 
the components arising from the bending of 1 for x — > 
— oo. Similarly to the AA cases, a phase slip (jump from 
branch to another) was usually found to be associated 
with a transition of the 1 texture. Textures with mirror- 



Figure |lj summarizes the critical mass currents J c for 
BB, AA, and AB junctions, plotted in the units defined 
in Appendix B. They have been calculated for several 
sizes in the range W/£gl,D /Cgl = . . . 10, which is 
certainly the most relevant one. For essentially larger 
dimensions the CPR's become strongly multivalued, and 
the aperture cannot be considered as a weak link. Fig- 
ure 15, on the other hand, is a phase diagram which re- 



lates the changeover between continuous and hysteretic 
CPR's, and some of the 1 texture transitions to well- 
defined regions of D, W plane. The results were calcu- 
lated using a lattice spacing of Ax/^gl — Ay/^GL = 0.5, 
and with a cutoff R c /£,gl = 20, but with no asymptotic 
corrections. In A A and AB, some amount of error will 
necessarily remain close to regions where 1 transitions oc- 
cur, since 1 is not allowed to vary freely beyond R c , which 
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FIG. 14: Critical currents for the four main phase combina- 
tions. The solid curves correspond to different widths W/£pr,. 
indicated by the numbering. BB is the same as in Ref. ESL 
but units differ by a factor of 10. The dashed lines show J c 
for apertures of fixed aspect ratio W : D = 2 : 1 — for BB 
two others are also shown. The asymptotic behavior at small 
sizes is the same for both AA cases. For AAjj, AAjj, and 
AB the J c 's become equal at large D. All J c 's are calculated 
for the most symmetric boundary conditions only. For other 
choices, the values (at small D) will tend to be smaller in BB 
and AA junctions, but larger in the AB junctions. The sharp 
features result from transitions of 1. They are explained in 
the text and in Fig. Ujil 



makes the texture somewhat too rigid. All of the results 
of Figs. ^ and ^| are for the most symmetric configura- 
tions only: in BB uj l = w , in AA d L — ±d R , and in 
AB — R^iSi. However, at least the narrow-aperture 
limit W/^gl % 3 is independent of any the boundary 
conditions. The regime with W/^gl 3 (for small D) 
is where the description due to Eq. (||) is supposed to be 
valid, whereas for larger apertures hysteresis begins to 
set in. This transition in behavior is seen in Figs. ^| and 

In Fig. |l4| the BB case is actually well-known. In Ref. 
it was calculated using a channel geometry, rather 
than the empty half-spaces on the two sides. The fact 
that the values in this reference are slightly larger may 
be just an indication of the restricted form of order pa- 
rameter which was assumed there. Thus we confirm the 
expectation that the critical currents should only depend 
on the properties of the weak link itself, not the way in 
which the current is driven through it. 

The other three panels of Fig. |lj show the interest- 
ing feature that if there is A phase on one of the sides, 
the critical currents will always approach the same val- 
ues for large D. They are thus determined by the inside 
of channel, and not the complicated 1 structures or the 
A-B interface at the ends of the junction. The transition 
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FIG. 15: Phase diagram for the 2D BB, A A, and AB junc- 
tions. For W/^gl < 3 the order parameter is strongly sup- 
pressed inside the hole, and the CPR's are continuous. For 
AAjj, this region is slightly larger. The letter 's' denotes re- 
gions where the A phase 1 field is symmetric, and 'a' regions 
where it is antisymmetric; 's/a' means that the configuration 
depends on Acj>. The sharp features in Fig. [[4] are always as- 
sociated with a change in the CPR branch which determines 
J c . In AAjj, the branch with 's' configuration, which exists 
for small Acj>, determines J c to the left of the solid line — the 
associated type of CPR is sketched in the inset. To the right 
of the line, J c becomes essentially independent of D, and the 
weak link approaches the "infinite slab" limit. Ji^] In AB, the 
regions marked by + or — denote whether J c arises from the 
(+) or the 7r/2 (— ) branch of the n periodic CPR's (cf. 
Fig. n2|) . To the right of the second transition marked by the 
solid line, the infinite slab limit is again reached. In AAjj no 
sharp branch transition exists, but the solid line has roughly 
the same interpretation. 



to the long-channel limit involves locking of 1 perpen- 
dicular to the channel walls (1 = ±y) at the middle of 
the junction, as in a parallel-plate geometry with plate 
separation VF.|32j For AAjj this happens smoothly, but 
in AAjj and AB there are rather sharp cusps in the J c 
curve beyond which J c is almost constant. In both cases 
this transition is due to a crossing over of the critical 
currents of two branches of CPR with different 1 configu- 
rations. For AA j j the configurations are the symmetric 
(e) and the asymmetric (d or f) ones in Fig. m, as dis- 
cussed in the previous section. Some of the associated 
CPR's were shown in Fig . [lC| , which corresponds to the 
W/£gl — 6 line of Fig. |14| — for small A(j) the sym- 
metric 1 configuration is the one to determine J c . The 
J c -transition lines are shown in the phase diagram, Fig. 
[la, although it should be kept in mind that their posi- 
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tions may depend on the way in which the numerics was 
done. For small D the behaviors of all four cases are quite 
different. For AAfJ. and AB J c tends to be suppressed, 
whereas for BB and AA|T it increases. In AB there are 
actually J c transitions at two different D values. Both 
are due to switching of J c between the (+) or tt/2 (— ) 
CPR branches, which are seen in Fig. [l^. The +/— re- 
gions are shown in Fig. [l5]. Notice that all these details 
are only valid for = dzR^Si. For other configurations 
the Fj un term is also finite, and the critical currents will 
in fact tend to be larger (for small D) than presented in 
the figures above. 

Also shown in Fig. |lj are the critical currents for 
some apertures of varying size, but fixed aspect ra- 
tio. These correspond essentially to Figs. [)] and |l2| for 
AA Tl and AB, and to equivalent (but not shown) fig- 
ures for BB and AA||. Small numerical differences are 
due to the different computational parameters. These 
curves are interesting because all the temperature depen- 
dence in the GL regime is in the length scale £,gl{T) ~ 
(1 - T/Tc)- 1 / 2 alone, and scaling D/£gl,W/£gl with 
factor / is thus equivalent to changing the temperature: 
I ~ (1 — T/Tc) 1 / 2 . The aperture in the Paris experiments 
reported in Ref. || was a slit of dimensions 0.18 x 2.6 fim 2 
in a 0.1 /im thick SiN membrane. The aspect ratio W : D 
for this is close to 2 : 1, which was therefore most fre- 
quently used in our calculations. The behavior at small 
I appears to follow a power law J c /Ja,b ~ l S - For all 
of the three BB cases the exponent is roughly the same, 
6 = 1.8. Similarly, for both A ATI and AA|| it is 5 = 2.2, 
whereas for AB 5 = 4.5. In BB the J c curves appear to 
be linear for large I, but this is only illusory. To point 
out another numerical accident, note that in the A ATT 
case the critical current transition appears to occur at 
exactly the aspect ratio W : D = 2 : 1. There seems to 
be no obvious reason for this. 

However, the striking similarity in the J c scaling be- 
havior of the AATI and AA|| configurations for small I 
appears very surprising at a first glance. Since in the pin- 
hole limit the critical current for AATI vanishes, should 
it not at least vanish faster as the smaller dimensions are 
approached? Two things should be realized here. The 
first is due to the 2D nature of our model: as the slit 
is made very narrow, eventually all but the z directional 
component of the orbital part of the order parameter 
are suppressed. The remaining coupling through a single 
complex component is direction-independent, and thus 
the behavior for both AA cases becomes similar. Second, 
even for a 3D aperture, the "pinhole" behavior would 
not be exactly reached in this limit, because the quasi- 
classical pinhole calculation relies on a nonlocal coupling 
through the "/ functions" . The present GL calculation 
is local, and assumes that the pairing potential is every- 
where nonzero inside the weak link — otherwise there 
could be no coupling at all. A different approach for a 
finite-size BB junction was recently adopted in Ref. |37j. 
These comments are also valid for the BB case, which 
explains the somewhat different result of Section 5 as 



compared with Ref. 



VII. CONCLUSIONS AND DISCUSSION 

In the sections above we have presented a p-wave 
Ginzburg-Landau analysis of a two-dimensional weak- 
link model. We have studied the current-phase relations 
and mapped the critical currents of the weak link, when 
it is placed between two bulk volumes of either superfluid 
3 He-A or 3 He-B. The BB case was studied before in sev- 
eral restricted calculations, [|J [| ^9) but here we were 
able to use much more general boundary conditions on 
the bulk phases. We showed that similar u tt states" as in 
BB pinholes (J, |j| exist also for small finite-size apertures 
in the GL regime. They are due to the second-order terms 
in Eq. @ and are thus associated with small critical cur- 
rents. On the other hand, we showed how this behavior 
gradually becomes more complicated in larger apertures, 
as the order parameter has more freedom to vary inside 
the hole. In this way, various kinds of 7r states may ex- 
ist together with either continuous or hysteretic CPR's 
which have relatively large critical currents. This is con- 
sistent with the experimental findings of Ref. where, 
apparently due to the uncontrollability of the bulk tex- 
tures, many kinds of CPR's were seen. As in Ref. § for 
the "anisotextural" effect in a large pinhole array, it was 
found that the hysteresis of CPR's tends to be smaller 
when the bulk boundary conditions are less symmetric. 
This enabled jumps to the it states which were unreach- 
able in symmetric cases. || Furthermore, we calculated 
the spin currents, and related their behavior to transi- 
tions between different order-parameter configurations in 
the aperture. We also studied AA and AB junctions, 
which have previously not been studied in any depth, at 
least not within the GL theory for large apertures. We 
analyzed their symmetries, and related transitions in the 
1 field to changes in forms of the CPR's and their criti- 
cal currents. In AB junctions, the usual phase-inversion 
symmetry was shown to be broken, which results in more 
complicated CPR's. In addition to all these, a hydro- 
static theory of the asymptotic phase corrections is pre- 
sented in Appendix A, and an analysis of the numerical 
relaxation method is given in Appendix B. 

We have also carried out some tests with the channel 
geometry used in Ref. and again we found the critical 
currents to match with those in Fig. 14, Also a geome- 



try with the weak link between two parallel flow channels 
was tested. It was found that the critical currents of a 
BB weak link are not affected by the parallel flow before 
the flow itself goes unstable by nucleating a vortex at 
some wall. However, these studies were not carried much 
further. What we have also not done above is an analy- 
sis of the effects of a strong magnetic field. The required 
magnetic GL terms would be relatively easy to include 
in the calculation, but they would add more dimensions 
to the already large parameter space. Actually there re- 
main even some untested phase combinations which may 
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be possible for junctions with small W/^ql > 3 and 
D > W: a BAB or ABA heterostructure. Such situ- 
ations may be obtainable at pressures close to Pq. This 
is because the A phase tends to be slightly more sta- 
ble in restricted "parallel plate" geometries, and on the 
other hand because there are some hysteresis effects re- 
lated to the equilibrium position of the A-B interface as 
the pressure is varied. |28| The BAB configuration was 
recently discussed as an explanation of some dissipation 
effects, |3S|] and this was, apparently, the initial reason for 
considering the BAAB structure in Rcf. ^. Although it 
should be possible to prepare such configurations in our 
program in some pressure regimes, their stability in a real 
experiment is questionable. Of course, the A phase could 
be stabilized more strongly inside the aperture by using 
a localized magnetic field. Also, since a self-consistent, 
general-temperature calculation of a finite-size aperture 
(with strong-coupling effects included) is still missing, 
no final conclusions can be drawn at this stage. Some of 
these issues may be worth further studies at some later 
time. 
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sumed to take the form: 

A(x) = fi s (x)i?A (0) (a;)[i? o (x)] T e iW+ ^ (x)) . (Al) 

Here <5</>(x) is a small phase correction and i? s (x) = 
R(S9 a (x)) a small spin rotation: R s v = 5^ — e^ va 59 a . 
In addition we have included a small orbital rotation 
i?°(x) = J?(5ej(x)). This is only needed in the anisotropic 
A phase to include corrections to the 1 texture, which 
are present even in the absence of currents. In deriv- 
ing the asymptotic corrections to 5<j), S9 a , Sei and the 
energy, we shall make the hydrodynamic approximation 
that A^ is everywhere in its bulk form. In the B phase 
this means that A^) — AbS^i and for the A phase 

A*jty = AA5 flx (S i y ± iS iz ). For convenience we shall de- 
note below the small-angle fields 5(f), S9 a , 5ei with </>, 9 a , 
6i, respectively. 

In this context we refer to all positions with respect to 
the origin (x,y) = (D/2,0). As a first approximation, 
all the small-angle corrections will be assumed to depend 
only on the radial distance r = \J x 1 + y 2 from this point, 
possibly with some simple dependence on the azimuthal 
angle (p = arctan(y/a;) also. 
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APPENDIX A: ASYMPTOTIC SOLUTIONS IN 

2D 

In this appendix we discuss the approximations used 
in the asymptotic regions between R c and i?oo of Fig. [I]. 
The need for these derives from a requirement to keep 
the numerical computation at a minimum, but at the 
same time have the bulk boundary conditions imposed 
as far away from the junction as possible. In 3D the 
boundary conditions (i?oo) could be taken all the way 
to infinity using this method, but in 2D some smaller 
value has to be chosen. Since in this case Roo is rather 
arbitrary, this procedure is to some extent only cosmetic, 
enabling a smoother transition of the order parameter 
at the R c cutoff. Thus, whenever this was not needed 
(like in evaluating the critical currents), the asymptotic 
corrections were not used. 

At Roc the boundary conditions are given by order 
parameters of the form A^ R — R^ R A^ LR (x)e 1<l ' I " R 
as shown above. However, since R and </) are con- 
stant, these states carry no mass or spin currents, apart 
from the spontaneous spin currents parallel to B-phase 
surfaces. pl[ The order parameters must be corrected in 
the asymptotic regions in order for them to conserve the 
currents coming from the weak link aperture. In this re- 
spect, the A and B phases can again be treated in the 
same way, and the asymptotics on each side can be as- 



1. Asymptotics in BW State 

We first deal with the case of B phase which is nu- 
merically simpler to handle than the A phase, because 
it is isotropic (meaning that it has an isotropic super- 
fluid density tensor). There is no need to worry about 
similar symmetry breaking effects as with the A phase 1 
vector (see below), and R° = 1. On the other hand, the 
"entanglement" of the spin and orbital parts of the order 
parameter makes their separation impossible, which com- 
plicates the analytic treatment. To facilitate the analysis 
of the gradient part of Eq. (||) , we first change to a spin 
basis where the constant R in Eq. (Al) is replaced by a 
unit matrix. This is done by noting that R S R = RR S , 
with R 3 ^ 



Sfw - e^ a 0' a , and 9' a = R^pOp, where the 



identity ei jk RjiRkm = Rin^nim is useful 

rent in this new basis is obtained with j^ la = R^pfp™ ■ 
However, to simplify notation we drop the primes below, 
remembering that they should appear on all quantities 
with Greek spin indices. With these definitions, the gra- 
dient energy (per unit length) can be written |p9| 



^iu-L. The spin cur- 



d 2 x [p s 



spin spin spin] 
' Paf3;ij V ai V /3j J 



and the mass and spin currents are 

j s = p s v s 

■ spin ^ spin spin 



(A2) 

(A3) 
(A4) 



where p s = (2m 3 /^) 2 2( 7 + 2)KA%, p s ^ = p s [(j + 
l)5ij8 al3 -(7- 2)5/3i5 a j -Saidpj], v s = (h/2m 3 )V(j) and 
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v a Pm = (fi/2wi3) V# Q . By making a variation in <j> and O 
and noting the symmetry P^"^ = p!, P rt"„- we see that 



Paf3y. 



the extrema of Fq are also found from the continu- 
ity equations V • j s = p s V ■ v s = {hp s /2m 3 )V 2 (j) = 
and V • jr 1 = = (W2ma)[( 7 + 1)V 2 0« 

— (7 — l)9o,9i0i] = 0. The latter set of equations can be 
directly diagonalized,[^2] but we assume the solutions to 
be radial right from the beginning: <fi(r, tp) = 4>{r) and 
6 a (r, ip) = 9 a (r) . Integration over angles in the energy 
[Eq. ( |A2| )] and minimization with respect to these inde- 
pendent phase fields yields, for example, the solution 



0(r, p) = <f>(r) = <p c 



(A5) 



These satisfy 4>(R C ) — 4> c and <j>(Roo) = 0. Similar ex- 
pressions are obtained for the spin phases 9 a (r, ip) . In the 
p s units defined in Appendix B, the ener gy corresponding 
to the solutions of the type of Eq. (A5) is 



F, 



2vr 



G 



4( 7 + 2) HR^/Rc) 



(7 + 2)(0 



c\2 



I7-I 

1 - -- — (i-s az 

2 7 + V 



and the asymptotic currents are 
J s = TT\n(R 00 /R c )- 1 (j) c 
jr n = ir(j + l)HRoo/Rc)- 1 



1 



17- 1 
27 + 1 



(1 - S az ) 



(A6) 

(A7) 
(A8) 



where no summation over a is implied. By fitting the 
numerically calculated currents J s , J^, pin (tra nsfo rmed 
with i? T to the "primed" spin basis) to Eqs. (A7) and 
( A8), one can solve the parameters (jf , 9 c a . Then one can 
use Eq. flA5| ) and the equivalent expression for 9 a in Eq. 
(Al) to update the the asymptotic order parameter, re- 
membering first to transform 9 c a back into the original 
"unprimed" spin coordinates. Th ere is an asymmetry in 
the spin parts of Eqs. (A6) and (A8): the z direction is 
slightly more rigid than the others. This results from the 
orbital-space (real-space) asymmetry of the 2D problem, 
which also affects the spin space in B phase. In the A 
phase this is not so, since there the spin and orbital parts 
of the order parameter are decoupled; see below. 



2. Asymptotics in ABM State 

In the A phase we assume the following conditions for 
the vector 1 = m x n: (i) 1 • w = 0, where the orbital 
rotation axis w is constant and in yz plane, so that (ii) 
w • z w 1, and (iii) 1 ■ x = cos e where the orbital rotation 
angle e is small. ^From condition (i) it follows that the 
phase <j> is well defined and acts as the superfluid velocity 



potential. Condition (ii) ensures that 1 is essentially in 
the xy plane with v s , and (iii) simply means that the 
variations from 1 = ±x are small. Using a linearized 
approximation for v s ^ n the A phase gradient energy fol- 
lowing from Eq. (^J) can then be written 



FA = 



1 



d 2 a 



[pijV si V s 



' Pap;ij u ai u j3j 



+ {h/2m 3 )v si C i:i (V x 1), + K S (V ■ I) 2 

+ K b \\x (Vxi)| 2 + i^(i-vxi) 2 ] 

and the mass and spin currents are 



(A9) 



Piji 



(fi/2m 3 )Cy(V x 1)3 



■ spin 

3 oii 



2m, 



spin spin 
~Paf3;ij V 0j J 



(A10) 

(All) 



where p^ = p±S i:j - (p± - p\\)klj, p± - p\\ = (7 - 1)/ (7 + 
1), p ± = (2m 3 /^) 2 2(7 + 1)KA% pf^ = 5 a0 p tJ v s = 
(h/2m 3 )V<f>, v^ in = (h/2m a )(d?d$> -6 a p)V9 , K s = 
K t - (^/2to 3 ) 2 ^/(7 + 1), K b = (^/2m 3 ) 2 p±7/(7 + 1) 



and C;, 



CiSi 



(C ± - Cu)kl r Here d c a 



the bulk spin vector. The first three terms of Eq. ( |A9| ) 
are just v s • j s /2 and v^ pin • j^ pin /2, but there arc three 
additional terms resulting from the 1 texture alone. Here 
the actual value of the C tensor will not be important 
to us, since with our assumptions the current component 
resulting from V x 1 is either divergenceless (1 in plane) or 
even vanishes (1 constant). As a result, we shall neglect 
the C terms also from the free energy, but keep all others 
initially. 

We now seek to minimize the energy [Eq. (|A9|)] with 
respect to the small-angle fields 0, 9 a and e. In principle 
these should be optimized simultaneously, but we shall 
separate the problems so that in optimizing <f> or 9 a we 
assume 1 = ±x, and in correcting 1 with e we assume 
that the velocities vanish. Assume, then, that we have 
1 = ±x, in which case the last three terms in Eq. ( |A9| ) 
vanish. By making a variation in <j> and 9 a and noting 
the symmetry j = P^p^j we see that the extrema 
of the first two terms of Fq are found simply from the 
continuity equations V j s = diPijV S j = and V • j^ pin = 
^Po^'ij^ft; 111 = 0- Now if we change to new coordinates 
x 1 = Xy/ (7 + l)/2, y' — y, then the superfluid density 
tensor p'^ will appear isotropic. In these coordinates the 
continuity equations are simply the Laplace equations 
V /2 (/> = V /2 #q, = 0. We attempt to solve these in polar 
coordinates r', ip' of the primed system, where the general 
solution is of the form 

<p{r',(p') = Aip' + Bhx(r'/a) 

OO 

+ 53(C„r' n + D n r'- n ) sin(V - a n ). (A12) 

n=l 

The solution we are looking for is determined by set- 
ting boundary conditions on the R c cutoff and at infinity. 
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FIG. 16: Schematic representation of a circulating current 
component (large arrows) in A phase due to the bending of 
the orbital anisotropy axis 1 (two-headed arrows). The small 
black arrows denote the symmetric, radially decaying current 
component. 



We assume that the currents coming from the junction 
flow radially as r' — > oo. This is described by the solu- 
tion proportional to lnr'. However, because the 1 tex- 
ture tends to align perpendicular to surfaces, it usually 
becomes nonsymmetrical close to the weak link (see Sec- 
tion 5). Since the superfmid density is larger perpendicu- 
lar to 1, most of the current will flow perpendicular to it. 
Therefore, the currents do not flow in and out of the junc- 
tion mirror-symmetrically with respect to the xz plane. 
Such a flow field can be thought to consist of two compo- 
nents: one which is symmetrical and radial, and another 
which circulates in front of the opening as illustrated in 
Fig. |l6|. For the circulating part of the mass current 
we test two boundary conditions: d r '4>(R c ) oc sin ip' and 
d r i<fi(R c ) oc sm2tp'. Both of these describe a flow which 
is outward above the x axis when it is inward below it 
and vice versa, but the distributions differ. The corre- 
sponding solutions are <f> = simp' /r' and <f> = sin 2^' /r' 2 
and the total velocity potential is thus of the form 



Hr'/Roo) 



,R 



+ c ' a ^sin(V), n = l,2. 



(A13) 

Here the coefficient <p c is again determined by equating 
the current obtained from Eq. (A10) by integration over 
R c (surface EDF in Fig. |l^) with the total current J s 
calculated numerically at the junction (surface AB). The 
coefficient (j> c ' a , on the other hand, is trickier. It could 
be determined by eq uatin g the current resulting from the 
second term of Eq. ( A13) over the upper part of the cut- 
off (surface DE) with the y-directional current J" cal- 
culated numerically over the surface CD. Note that in 
the case n = 1 the circulating parts of the currents over 
surfaces DE, Doo and DF are equal, and hence no net 
y directional current is introduced. This is not so for 
n = 2, and, in fact, in this case the phase field gives an 
unphysical-looking phase gradient perpendicular to the 
wall. However, our approximation does not take into ac- 



count the fact that the order parameter is suppressed at 
the wall and (at least for a diffusive surface) the super- 
fluid density vanishes there. Therefore no cu rrent should 
flow through it. The phase correction of Eq. ( A13 ) should 
only be seen as a variational ansatz, and in practice we 
find that n — 2 gives a better overall fit with lower total 
energies than n = 1. Similar expressions and conclusions 
exist for the spin velocity potentials 6 a . 

We now give the results for the current vs. phase cor- 
rection relationships. Let us remind that in deriving 
these we have assumed that 1 = ±x in the asymptotic 
region. In the p± units defined in Appendix B, the gra- 
dient energy of the asymptotic solution is 



2 V 7 + 1 1 HRooARc) 



E 



ln(iW#c) 



2n 

^ 2n 



(R c /RocY n Wa a ) 2 



while the currents are 



Js 


= V2/(7H 




rspin 


= 7r v /2/( 7 H 


-lJMiJoo/iJc)- 1 ^ 




= 71^2/(7 H 


- l)0 c ' a 



J s /^ a = ^2/(7+1)^ 



(A14) 

(A15) 
(A16) 
(A17) 
(A18) 



Again, by fitting the numerically calculated currents to 
these expressions, one can solve the parameters <fi c , 8 c a . 
The above methods of estimating the "asymmetric" cur- 
rents are very unreliable, and in practice the corrections 
(j> c,a and 0^ a were not used. They are presented here for 
completeness. 

In addition there is the small correction to the direc- 
tion of 1, which is mostly due to the complicated surface 
structure. There is also a small tendency for the 1 to 
align with ±v s at high superfiow velocities, but here we 
assume v s = v^ pm = 0. We do not require 1 to vary only 
in xy plane, as long as it is rotates only around the fixed 
axis w which is in the yz plane. The minimization of 
the last two terms of Eq. (A9) is done by assuming that 



l(r', ip') = ±i?(w, e(r', ip'))*., so that 1 = ±x for r' — > oo 
but also at the walls (ip' = ±7r/2). This can be described 
with the variational ansatz e(r',i?') = e(r') cos ip'. If we 
further assume that (K Sjt + K b )/2 < K b (^ + l)/2 in 
order to drop some terms with y and z derivatives, the 
energy density resulting from the bending of 1 is just 
2f l g = Kb[(d x e y ) 2 + (d x e z ) 2 }, where e = ew. Minimiza- 
tion of this with the above ansatz gives e(r') e c (R c /r'), 
and the energy 



37T7 

~16 



7 + 1 



[l-(#c/i?oo) 2 ](e c ) 2 



(A19) 



Here (e c ) 2 = (eg) 2 + (el) 2 , and » l^x' = R c ,y' = 
0), These cannot be obtained from any conservation law, 
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but instead we extract them from the numerical solution 
by requiring continuity of l y and l z and their x deriva- 
tives along the x axis. The 1 vector can be obtained 
numerically from the A matrix inside the R c region with 
k = tijk Im(A* iJ 4 / _ t fc)/2, as shown in Fig. |[ 

APPENDIX B: NUMERICAL METHODS 

This appendix considers in some detail the implemen- 
tation of our numerical relaxation methods. We skip the 
easier part of calculating the ID order parameter (x) 
for a planar wall, since it is only needed as a boundary 
condition for the 2D calculation. Of course, the same 
techniques can be used for that also, but in practice the 
number of variables in the ID minimization is so small 
that perfecting the method is not of vital importance. At 
all solid surfaces, we only consider the diffuse-scattering 
boundary condition A^ = 0, see Section 3. Below we 
assume that the boundary conditions on the minimiza- 
tion region are fixed, and comment on the addition of the 
asymptotic corrections (Appendix A) in the last subsec- 
tion. For the purpose of numerics, we must consider all 
quantities to be dimensionless, and we first discuss the 
unit reductions. 



1. Units 

A natural way to scale any quantity L with the units 
of length is L = L/^gl^ and a dimensionless order pa- 
rameter is obtained with A = A/ Aa.b, where Aa,b are 
the bulk gaps of A and B phases, respectively. In the A 
phase A 2 A = a/[2(2/3 245 )], where /3 245 = P2 + P± + A>, 
and for the B phase A 2 B = a/[2(3/?i2 + /?345)], where 
P12 = Pi + P2, #345 = 03 + Pi + As- Physically motivated 
units for the energies and currents are based on the the 
bulk superfluid densities p^ = (2m 3 /K) 2 2(^f + 1)KA 2 A 
and p B = (2m 3 /h) 2 2(j + 2)KA' B (see appendix A). Using 
these we define the following units for the currents and 
energy (per unit length). In A phase J a — (h/2m 3 )p±, 
J s / n = {h/2mz) 2 p ± , and E A = {h/2m 3 ) 2 p ± . For 
B phase Jb — (fr/2m 3 )p s , J B pm = (h/2m 3 ) 2 p s , and 
Eb = (h/2m 3 ) 2 p s . It is these units which are used in 



all the figures of this paper, and sometimes we denote 
Js = Js/Ja,b, JT n = J spin /J s /B\ and F = F/E a ,b- 
Note that for weak coupling (see Section 3), the A and 
B phase superfluid densities p s and p± are equal, and so 
are the units. For cases where both A and B phase are 
involved simultaneously, we choose to use the B phase 
units. 

However, for numerics, a more convenient reduction of 
the energy is given by / = / /[aA AB /2], and then the 

reduced GL parameters are given by p i = (2A A g /a)fti, 
a = 2 and K — 2. Using the Aa,b we thus have in the A 
phase P i = A/(2/3 2 45) and in the B phase P i = Pi/ (3/3i2+ 
/?245). In the formulas below, these unit reductions are 



assumed to have been performed, but we do not write 
any tildes or overlines on the symbols. 



2. The Minimization Problem 

The 2D numerical scheme is based on an optimiza- 
tion of the free energy Fq [Eq. Wj] or, which is more-or- 
less equivalent, solution of the corresponding nonlinear 
Euler-Lagrange field equations. In fact, the method to 
be discussed can be applied (at least approximately) for 
solving a more general class of nonlinear equations which 
have the form G /i j(A(x)) = 0, where p, i — x, y, z. Nev- 
ertheless, we shall consider the the problem as that of 
optimizing (minimizing) the GL energy functional Fq [A] 
in the region £1 of Fig. [I]. More mathematically speak- 
ing, for the gi ven fixed boundary conditions A L ' R [of the 
form of Eq. (|Il|) 1 on T L > R , we wish to find F j[A L , A R ] 
= min^ Fq [A] . As discussed in Appendix A, the full 
numerical minimization is done only inside the region 
bounded by the R c cutoffs, but the computational region 
can be effectively increased by making use of the asymp- 
totic corrections between R c and i?^. 

In the following, all entities A with the 3x3 matrix 
structure A^i where p,i = x,y,z will be referred to as 



spin-orbit matrices, or spin-orbit fields. For simplicity of 
notation, let us define the following rotational invariants 
which are functions of the spin-orbit matrices A, B,C, D: 



fa(AB) 


= ReTr(AB T *), 


(Bl) 


fi(A,B,C,D) 


= RcTr(A*B T *) Tr(CL> T ), 


(B2) 


f 2 (A,B,C,D) 


= RcTr(AB T *) Tr(C*L» T *), 


(B3) 


f 3 (A,B,C,D) 


= Rc Tt(AB t C*D' t * ) , 


(B4) 


U(A,B,C,D) 


= ReTr(AB T *C£> T *), 


(B5) 


f 5 (A,B,C,D) 


= Rc Tt(AB t *C*D t ) , 


(B6) 



g 1 (A,B)=Red i A; i d j B lij , 
g 2 (A 7 B)=Red l A; j d l B M , 
g 3 (A,B)=Red i A* tlj d j B fli . 



(B7) 



The invariants have several symmetries, some of which 
are listed below: 

f a (A,B)=f a (B,A)=f a (A*,B*), (B8) 

h a {A, B, C, D) = f lt2 (A, B, D, CO, (B9) 

/ 3 ,4, 5 (A B, C, D) = f 3A5 (B, A, D, C), (BIO) 

/i_ B (4 B, C, D) = A_ 5 (C, D, A, B), (Bll) 

f!- 6 (A*,B*,C, D) = A_ 5 (A B, C*,D*), (B12) 

01,2,3(4 B) = 91 , 2 , 3 {B, A) = g 1A3 {A*,B*), (B13) 



It is worth noting here that symmetry of Eq. ( |B13| ) of the 
gradient terms is not necessarily valid in a discretized 
form, unless symmetric difference approximations are 
used for all derivatives. Using the definitions of Eqs. 
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( piyB7| ), the GL energy functional [Eqs. (||) and (§)] can 
now be rewritten as follows 

r 3 

F[A,A*]= / d 3 x[- af a (A,A) +Y,K tgi (A,A) 

J z=l 
5 

+ ^A/i(AAA^)]- (B14) 

i=l 

For the purpose of numerics, it is very helpful to con- 
sider F as a functional of both A and A* independently, 
as shown he re explicitly. Note that the third gradient 
term in (B7) can be integrated by parts and added to 
the first term, if a separate surface energy term is intro- 
duced. Then, if the surface and the associated boundary 
conditions are symmetric enough, the surface term may 
vanish altogether, thereby reducing the numerical effort 
somewhat. Since this is not always the case, we shall not 
assume it anywhere. 



3. Line Searching in Minimization 

The basic component of any efficient minimization al- 
gorithm is a method for doing line searches, i.e., doing 
one-dimensional minimizations of the function in some 
given direction. More precisely, one wishes to minimize 



h(X) = F[A + XD,A* + XD*} 



(B15) 



with respect to the real parameter A, given the "starting 
point" A and a "search direction" D. The extrema of 
h(X) are found from the zeros of its derivatives, i.e., by 
solving 



h'{X) = 2 Re (G(A + XD),D) = 0, 



(B16) 



where G(A(x)) = SF[A, A*]/5A*(x), and we have de- 
fined a scalar product of two spin-orbit fields A and B as 
follows: 

(A,B) = J d 3 a;Tr(A(x)B T *(x)). (B17) 

It is easy to check that this is indeed bilinear, satisfies 
(B, A) = (A, B)* , and that (A, A) > 0, with (A, A) = 
only for A = 0. The scalar product also gives a natu- 
ral definition of a norm for the spin-orbit field: \\A\\ — 
y (A, A). Note that, with respect to the matrix indices, 
this is the so-called Frobenius norm. 

The interesting point is that, since the GL free energy 
is only of fourth order (quartic) in the order parameter, 
the function h(X) is simply a fourth order polynomial 
with real coefficients 



h(X) = aA 4 + 6A 3 + cA 2 + dX + e. 



(B18) 



The fourth-order coefficient is clearly always positive, so 
that the graph of the function must look like the curve in 
Fig. [It]. The number of minima is 1 or 2, but only those 



h(X) 




FIG. 17: A fourth order polynomial (solid line) and a local 
quadratic approximation to it (dashed line) around A = 
(dotted line). 



with A > will be of interest here. If one is now willing to 
take the small trouble of calculating the coefficients a to 
e, then the exact minima of h(X) are immediately found 
by solving the real roots of a third-order equation h'(X) = 
0. For the conjugate gradient method the accuracy of line 
searches is important, at least in theory, and we list below 
the formulas for a to e for the GL free energy discussed 
above: 

5 

a = J d 3 x[Y,Pifi(D,D,D,D)] (B19) 

i—1 



5 

b= I d 3 x[4j2Pifi(D,D,D,A)} (B20) 

i=l 

5 

= / d s x{2j2 Wi(D, D, A, A) +fi(D, A, D, A) 

J i=l 

3 

+ /,(£>, A, A, D)] - af a (D, D) + ^ K l9l (D, D)} 

(B21) 

5 

d= d 3 x{4 : Y / foMD,A,A,A)-2af a (D,A) 

i—1 

3 

+ Y / K>{ 9 i(D,A)+g l (A,D)}} (B22) 



e = F[A,A*} 



(B23) 



Here the last term of d can now be simplified with Eq. 
(B13), but only with the proviso that the difference ap- 
proximations are symmetric. The coefficient e (free en- 
ergy at position A) is just a constant and unnecessary for 
the line searching procedure. 
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The numerical calculation can be further simplified by 
assuming the coefficients a and b to be zero. This is 
equivalent to a local quadratic approximation of the en- 
ergy functional, which is reasonable if the center point 
A is already clo se to the minimum of F (see Fig. 17). 
Note that Eqs. (B19-B22) for the c oeffic ients a to d can 
also be obtained directly from Eq. (B16), given the GL 
equations G(A) — 0. This is indeed the only way, if 
no such functional exists whose Euler-Lagrange equation 
this G(A) = would be. Of course, in that case there is 
no guarantee that the CG iteration will converge at all. 



4. The Conjugate Gradient Method 

Perhaps the simplest approach for the minimization 
of a nonlinear function(al) is the steepest descent (SD) 
method. In this method, on each (fcth) step of iteration, 
the line search direction is simply chosen to be that of 



the negative gradient: D k 



-G k , where G k 



G(A k 



The line search condition [Eq. (B16)] then implies that 
(Dfc, -Dfc-i) = 0, i.e., successive search directions are per- 
pendicular to each o ther with respect to the appropriate 
scalar product [Eq. (B17)]. This leads to the well-known 
zigzag-type trajectory, which in many cases results in ex- 
tremely slow convergence. p0| The situation can usually 
be helped somewhat by skipping the costly line search 
procedure altogether, by choosing A k+1 — A k — c ■ Gk, 
where c is a given small positive constant. This proce- 
dure is in fact equivalent to a simple ("forward Euler") 
time-discretization of the time-dependent GL (or TDGL) 
equations, where c is proportional to the time step and 
the strength of dissipation. [^5| By following the progress 
of this iteration one may thus get a rough idea of the 
dynamics of the order parameter when, for example, a 
phase slip in a weak link occurs. 

At least in the case of the GL equations, a significant 
speedup of convergence can be achieved by taking advan- 
tage of the knowledge of previous search directions. In 
the standard conjugate gradient (CG) methods fhl the 
iteration step fc is taken as follows: 

Conjugate gradient iteration step 

• D k = -Gk+{ik-D k - l 

• Do a line search of h(X) = F[A k + XD k ,A* k + 
XD k ] to find its (smallest) positive minimum 
point A 

• Set A k+1 = A k + X-D k 

The constant (3 k can be calculated in one of several ways. 
Two of the most commonly used are the Fletcher-Reeves 
choice (3 k = (G k ,G k )/(G k -x,G k ^i), and the Polak- 
Ribiere choice f3 k = (G k - G k - 1 ,G k )/(G k - 1 ,G k - 1 ). 
The latter can be slightly more efficient, but it has the 
memory-economical disadvantage that the gradient G k _\ 
of the previous round must also be stored. The method 
is usually started with Dq = —Go which is guaranteed 



to be a direction of descent. During the iteration the 
process can also be restarted every once in a while. This 
is because it is possible that inefficient search directions 
will begin to be generated after a few rounds of mini- 
mization, if the function(al) deviates very strongly from 
a quadratic form. Restarting should be done at least if 
the line search fails to find a positive A, and possibly on 
every ATh step or so, where N is the total number of real 
variables in the discretized order-parameter field. |^0| 

It seems that in practice neither the accuracy of the 
line search, nor the choice for the f3 k s is very important, 
at least if the method is sometimes restarted. The most 
important part is really the basic philosophy of using the 
previous search direction as described above. 

One of the benefits of exact line searching is that 
it provides the small change in energy between steps 
very naturally through AF k — aX 4 + bX 3 + cX 2 + dX, 
so that e, the energy itself, does not need to be cal- 
culated on every round. This appears to be a more 
accurate way to obtain the change than application of 
AF k = F[A k ,A* k ] - F[A k _x,A* k _^ after completion of 
each iteration step. This is an important point when the 
limit of available floating point accuracy is eventually ap- 
proached: one should always have AF k < for a descent 
direction D k , and it is wise to check that this is indeed 
the case. However, a mistake would be done by assuming 
that F[A n , A* n ] = Y,tZo AF k + F[A , A*}, since this will 
involve summation of numbers whose orders of magni- 
tude are very different from each other. Sufficiently close 
to the minimum an F calculated in this way may no 
longer change at all during the iteration, even if the the 
gradient is still nonzero. The change in energy should 
therefore be used as a judge of convergence only with 
care. 



5. Comparison of Methods 

Here we compare briefly the convergence of the dif- 
ferent methods. We use an example calculation in the 
geometry of Fig. 0, with R C /£ GL = 25, W/fe L = 10, 
D/£gl = 6, with no asymptotic corrections (i?^ = R c ), 
and B phase with Cj l = lj r on both sides. The dis- 
cretization lattice spacing was chosen to be Ax/^gl — 
Ay/iGL — 1- Figure [l^ shows plots of the length of the 
gradient on the fcth iteration round, ||Gfe|| = ||G(j4fc)||, 
versus the round number fc. Both the Fletcher- Reeves 
and Polak-Ribiere variants of the CG method are com- 
pared to the SD method. The superiority of either form 
of CG is evident from this figure: while the CG methods 
reach the level ||G fc || w 10~ 4 for k = 200 . . .300, the SD 
method requires almost fc = 4000 for it. The real time 
needed by SD is thus more than tenfold compared to CG. 
This does not reveal anything rigorous about the conver- 
gence classes, but for the CG methods the convergence 
must be superlinear, and probably close to quadratic. |40| 
The scaling properties (behavior of convergence with in- 
creasing problem size) were also not tested. As an inter- 
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arbitrary way to form a vector. In the case of a compu- 
tational region looking like the one in Fig. [j] and not a 
rectangle, there is no obvious way to do such ordering. 
Any choice would require one to program obscure and 
lengthy routines for transformations between the "xy" 
and vector orderings. The CG method requires no such 
thing, and all field variables can be iterated in any order 
on each iteration step, independently of each other, (iv) 
The explicit use of A and A* as the independent vari- 
ables instead of Re A and ImA is also simplest in CG 
type methods. The practical necessity of this is clear 
from the complicated structure of the free energy func- 
tional. Points (iii) and (iv) are also why an application of 
most ready-made library routines is not very attractive 
for problems in general p-wave Ginzburg-Landau theory. 



FIG. 18: Gradient length versus number of finished iteration 
rounds on a log- log scale. The solid and dashed lines are for 
the Polak-Ribiere and Fletcher-Reeves variants of the conju- 
gate gradient method, respectively, and the dash-dotted line 
for the steepest descent method. All methods use exact line 
searches, so that the real time spent on one round is practi- 
cally the same for all of them. For a tolerance level of 10~ 4 — 
where the steepest descent method was terminated by hand 
— greater than order-of-magnitude differences in numbers of 
required iteration steps are visible. The conjugate gradient 
methods were restarted on round 200, which did not affect 
the convergence in any noticeable way. 



esting detail, it was noticed that the A phase tended to 
converge faster than B phase. 

Another possibility for obtaining faster convergence 
compared to the SD algorithm would be to use quasi- 
Newton (secant) methods, which utilize knowledge of the 
Hessian matrix or some approximation to it.pO| How- 
ever, in our example problem the CG method has several 
advantages over these methods also: (i) While analytic 
calculation of the second derivatives of F is possible, the 
Hessian matrix can easily become huge and take up a lot 
of memory. This is not so severe in approximate schemes 
(such as BFGS) which do not store the whole Hessian 
matrix, (ii) The Newton method has an inconvenient 
property of being unstable, if the initial guess is not cho- 
sen close enough to the minimum. To avoid problems, 
special stability checks are needed, (iii) Most inconve- 
niently of all, the construction of the Hessian matrix and 
solving the related linear system of equations would re- 
quire the programmer to order the field variables in some 



6. Updating the Asymptotic Values 

Finally we consider the task of fitting the asymptotic 
solutions of Appendix A to the numerical solutions ob- 
tained by the above methods for A L,R fixed on r = R c . 
This could be done in several ways, but here we choose a 
very simple one, and describe it only in therms of the B 
phase. The basic idea is iterative. Start with cjf = d c a = 
in Eq. (A5), an thus no phase c orre ctions in the the 
asymptotic order parameter [Eq. (|Al|)]. Then (i) min- 
imize F inside the R c cutoffs, (ii) calculate the currents 
J s and J^ pin in the junction, (iii) insert them in E qs. ( |A7| ) 
and ( [AS) ) to get new </> c and 6^, (iv) update Eq. ( |Al"| ) ac- 
cordingly, and then repeat from (i) until satisfactory con- 
servation of the currents over the R c cutoff is achieved. If 
we denote by "out" the current in the asymptotic region, 
and AJ S = J° ut — J s then the requirement was usually 
\AJ S \ < 10 -4 , and similarly for the spin currents. 

Note that A J s cx dFn/d(/) c , where Fq is the energy in 
the whole for given </> c , 6^j, and thus the process (i-iv) 
is seeking the minimum of Fq . Since this is obviously a 
very inefficient method, the benefit of fast convergence 
by the conjugate gradient method in the r < R c part of 
£1 is lost if this iteration is continued very far. Fortu- 
nately, the reason for taking the corrections into account 
is mostly cosmetic in practice, and only a few (less than 
10) iteration rounds may already lead to a reasonably 
continuous order parameter and current across the R c 
cutoff. In the case of A phase, also the correction of 1, 
and possibly the asymmetric phase corrections, should 
be updated on each round. See Appendix A for details. 
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